In this document, I want to provide a write-up of how this model is similar to and different from the versions previously implemented. This model is meant to provide assembly of multiple sites at the same time, that may, or may not, be connected in some fashion. Along the way, I will mention some of the input and output format that I am expect as documentation. At the end, I will be presenting some preliminary results. I especially want to use this as a vehicle to think about how to analyse those results.
library(dplyr) # Data Manipulation
Warning: package ‘dplyr’ was built under R version 4.1.1
Attaching package: ‘dplyr’
The following objects are masked from ‘package:stats’:
filter, lag
The following objects are masked from ‘package:base’:
intersect, setdiff, setequal, union
library(ggplot2) # 2-D Plot
Warning: package ‘ggplot2’ was built under R version 4.1.1
library(plotly) # 3-D Plot
Warning: package ‘plotly’ was built under R version 4.1.1
Registered S3 method overwritten by 'data.table':
method from
print.data.table
Registered S3 methods overwritten by 'htmltools':
method from
print.html tools:rstudio
print.shiny.tag tools:rstudio
print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
method from
print.htmlwidget tools:rstudio
Attaching package: ‘plotly’
The following object is masked from ‘package:ggplot2’:
last_plot
The following object is masked from ‘package:stats’:
filter
The following object is masked from ‘package:graphics’:
layout
library(ggfortify) # used for biplots of PCAs
Warning: package ‘ggfortify’ was built under R version 4.1.1
library(vegan) # Ecological analysis mega-package
Warning: package ‘vegan’ was built under R version 4.1.1
Loading required package: permute
Warning: package ‘permute’ was built under R version 4.1.1
Loading required package: lattice
This is vegan 2.5-7
library(RMTRCode2)
# https://stackoverflow.com/a/7172832
ifrm <- function(obj, env = globalenv()) {
obj <- deparse(substitute(obj))
if(exists(obj, envir = env)) {
rm(list = obj, envir = env)
}
}
First, we load up the preliminary results. In this case, we are loading a system in which 34 basal species and 66 consumer species form the pool for 10 unconnected environments. The pool and interaction matrix were assembled with the default parameters from Law and Morton’s 1996 work.
load(file.path(
"..", "experiments", "MNA-FirstAttempt-Result-Env10-None.RData")
)
toRemove <- result$Abundance <= result$Parameters$EliminationThreshold
result$Abundance[toRemove] <- 0
ifrm(toRemove)
In total, 9320 events were used in these environments, with the species and environment invasion both randomly assigned. The number of arrival and extinction events were controlled to both be half of this number. We chose this number due to the coupon collecting problem. In particular, we use the result that the probability of encountering each species is bounded: \[\text{Pr}(\text{Draws} < n \log_{e} n + c n) \rightarrow \exp(-\exp(-c)) \text{ as } n \rightarrow \infty\] where \(n\) is the number of species and \(c\) is a constant. For our purposes, we choose \(c = 5\) so that we have a probability of about \(99.3\%\) of seeing each species in each environment. In practice, we failed to observe 2 species-environment combinations. Notably, nearly every species had at least one successful invasion; 17 did not.
The initial abundance was set to be 4000 times the elimination threshold, in line with work on minimum viable populations (Traill et al. 2007). The elimination threshold is admittedly more arbitrary, since it sets an effective individual-area relationship. For this calculation, I set it to \(10^{-4}\), in line with our previous calculations. This is large enough to avoid numerical difficulties from precision, while being low enough to represent a decent sized region.
Since each population is assembling simultaneously, I chose to use exponential waiting times for the inter-species arrival and extinction times. Note that these rates are shared between species and environments, but arrival and extinction are fully independent of each other. Species and environment affected in each event was chosen uniformly at random. The question then is how to set the rate.
To set the rate in this case, I chose to set it to the largest eigenvalue magnitude of the per-environment interaction matrices. This magnitude corresponds to the strongest response we can see from the interaction matrix and, if the interaction matrix is a good approximation for the Jacobian around a stable fixed point (which is not guaranteed), indicates the characteristic time scale of the decay to equilibrium. Hence, (overall) arrival rates and (overall) extinction rates should happen on the same timescale as the (largest) dynamics in the system. Since there are 10 environments, we should then expect that 10 characteristic time scales, on average, should occur in between arrival events in the same environment.
With 10 environments, it is probably not helpful to check 10 individual abundance curves, but looking at the first one might be helpful.
LawMorton1996_PlotAbundance(result$Abundance[seq(from = 1,
to = nrow(result$Abundance),
by = 10), c(1, 2:101)]) -> obj;
obj + ggplot2::scale_y_log10() + ggplot2::guides(color = FALSE)
Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
Warning: Transformation introduced infinite values in continuous y-axis
Every vertical line is a species introduction or extinction by neutral dynamics.
Perhaps more intriguing might be some sense of the biodiversity that we have in each system. We break the abundance results up by environment, then calculate the number of non-zero abundance curves at each time point. We also calculate the Shannon entropy (reminder: higher entropy means more uncertainty which means a flatter distribution).
Diversity <- lapply(
1:result$NumEnvironments,
function(i, abund, numSpecies) {
time <- abund[, 1]
env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
richness <- rowSums(env != 0)
abundSum <- rowSums(env)
entropy <- env / abundSum
entropy <- - apply(
entropy, MARGIN = 1,
FUN = function(x) {
sum(ifelse(x != 0, x * log(x), 0))
})
species <- apply(
env, MARGIN = 1,
FUN = function(x) {
toString(which(x > 0))
}
)
data.frame(Time = time,
Richness = richness,
Entropy = entropy,
Species = species,
Environment = i)
},
abund = result$Abundance,
numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)
Diversity <- dplyr::bind_rows(Diversity)
Diversity <- Diversity %>% dplyr::mutate(
Evenness = Entropy / log(Richness)
)
ggplot2::ggplot(
Diversity,
ggplot2::aes(
x = Time,
y = Richness,
color = factor(Environment),
alpha = ifelse(Environment == 1, 1, 0.3)
)
) + ggplot2::geom_line(
) + ggplot2::geom_line(
data = Diversity %>% dplyr::group_by(
Time
) %>% dplyr::summarise(
Richness = mean(Richness)
),
mapping = ggplot2::aes(
x = Time,
y = Richness
),
color = "black",
inherit.aes = FALSE
) + ggplot2::guides(
alpha = "none"
) + ggplot2::scale_color_discrete(
"Environment"
)
So richness hovers around a similar regime throughout the majority of the simulation. Note in this plot that we have emphasised one environmental curve and superimposed the mean in black. We do manage to reach heights of 11 species in one environment, but these heights are shortlived. Instead, we seem to observe a (time and environment averaged) value:
mean(Diversity$Richness)
[1] 4.458893
(If we consider the first 10,000 time units as burn-in, we instead see a value of
mean(Diversity$Richness[Diversity$Time > 10000])
[1] 4.587217
ggplot2::ggplot(
Diversity,
ggplot2::aes(
x = Time,
y = Entropy,
color = factor(Environment),
alpha = ifelse(Environment == 1, 1, 0.3)
)
) + ggplot2::geom_line(
) + ggplot2::geom_line(
data = Diversity %>% dplyr::group_by(
Time
) %>% dplyr::summarise(
Entropy = mean(Entropy)
),
mapping = ggplot2::aes(
x = Time,
y = Entropy
),
color = "black",
inherit.aes = FALSE
) + ggplot2::guides(
alpha = "none"
) + ggplot2::scale_color_discrete(
"Environment"
)
Warning: Removed 8043 row(s) containing missing values (geom_path).
Warning: Removed 2297 row(s) containing missing values (geom_path).
Entropy has a similar, but highly erratic, behaviour. If one tries to follow one of the entropy curves, then one sees that they have fairly substantial periods of almost smooth behaviour followed by suddenly very noisy behaviour, and noise seems to be the dominant mode if one tries to examine the low alpha environment in the background. There are some easy to make predictions. Extinctions reduce the entropy in the system, as you become more certain about what remains. Analogously, arrivals increase the entropy. We can probably better see the relationship beyond these principles by plotting entropy against richness and connecting observations by environment and time.
# ggplot2::ggplot(
# Diversity %>% dplyr::filter(Environment < 4),
# ggplot2::aes(
# x = Richness,
# y = Entropy,
# group = Environment,
# color = Time
# )
# ) + ggplot2::geom_path(
# ) + ggplot2::guides(
# alpha = "none"
# )
plotly::plot_ly(data = Diversity %>% dplyr::filter(Environment < 2),
x = ~Richness, y = ~Entropy, z = ~Time, type = "scatter3d",
mode = "lines", opacity = 1, line = list(color = ~Time))
It seems quite hard to tell, but there does not appear to be any particular orientation (clockwise, counter-clockwise) or similar pattern here.
ggplot2::ggplot(
Diversity,
ggplot2::aes(
x = Time,
y = Evenness,
color = factor(Environment),
alpha = ifelse(Environment == 1, 1, 0.3)
)
) + ggplot2::geom_line(
) + ggplot2::geom_line(
data = Diversity %>% dplyr::group_by(
Time
) %>% dplyr::summarise(
Evenness = mean(Evenness)
),
mapping = ggplot2::aes(
x = Time,
y = Evenness
),
color = "black",
inherit.aes = FALSE
) + ggplot2::guides(
alpha = "none"
) + ggplot2::scale_color_discrete(
"Environment"
)
Warning: Removed 10986 row(s) containing missing values (geom_path).
Warning: Removed 2831 row(s) containing missing values (geom_path).
Evenness helps highlight that this is a high variance process but with a relatively constrained mean.
We can, of course, flip the idea on its head. Instead of examining the diversity of species within environments, we can look at the diversity of environments occupied by species. Since very few species end up occupying environments, we just look at richness. Unfortunately, this is quite a memory exhaustive task.
EnvDiversity <- lapply(
1:((ncol(result$Abundance) - 1) / result$NumEnvironments),
function(i, abund, numSpecies) {
time <- abund[, 1]
env <- abund[, 1 + i + numSpecies * (1:result$NumEnvironments - 1)]
richness <- rowSums(env != 0)
abundSum <- rowSums(env)
environments <- apply(
env, MARGIN = 1,
FUN = function(x) {
toString(which(x > 0))
}
)
data.frame(Time = time,
Richness = richness,
Abundance = abundSum,
Species = i,
Environments = environments)
},
abund = result$Abundance,
numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)
EnvDiversity <- dplyr::bind_rows(EnvDiversity)
ggplot2::ggplot(
EnvDiversity %>% dplyr::filter(Richness > 1),
aes(x = Time, y = Richness, color = Species)
) + geom_point(
alpha = 0.01, size = 3
) + guides(
color = "none"
)
One immediately interesting trend here is that very few species are present across more than 5 environments at a given time. Indeed, only
unique(EnvDiversity$Species[EnvDiversity$Richness > 5])
[1] 12 14 22 30
are ever present in more than 5 environments at once. We can also examine how long these periods occur for by species by tabulation. We block times so that entries are all of the same unit length, and round the average richness during the given time unit.
with(EnvDiversity %>% mutate(
Time = floor(Time)
) %>% group_by(
Time, Species
) %>% summarise(
Richness = round(mean(Richness)),
.groups = "drop"
),
table(Species, Richness))
Richness
Species 0 1 2 3 4 5 6 7
1 65970 14384 439 0 0 0 0 0
2 34776 31378 11352 3287 0 0 0 0
3 63206 13000 1767 2820 0 0 0 0
4 46367 21078 11443 1617 288 0 0 0
5 30648 31553 16542 1341 709 0 0 0
6 72061 8732 0 0 0 0 0 0
7 75522 5271 0 0 0 0 0 0
8 64759 15508 526 0 0 0 0 0
9 77919 2874 0 0 0 0 0 0
10 17817 38679 20192 4105 0 0 0 0
11 52666 24892 3235 0 0 0 0 0
12 9669 2676 9781 26641 18597 9922 2294 1213
13 64640 15703 450 0 0 0 0 0
14 12094 21314 14556 16445 15148 931 305 0
15 68564 10732 1497 0 0 0 0 0
16 69690 10903 200 0 0 0 0 0
17 30636 33714 15505 938 0 0 0 0
18 69214 11322 257 0 0 0 0 0
19 62460 17599 734 0 0 0 0 0
20 74867 3685 1849 392 0 0 0 0
21 47906 25702 4185 3000 0 0 0 0
22 3180 13477 20151 21766 16740 5305 174 0
23 65137 14940 716 0 0 0 0 0
24 65792 12313 2688 0 0 0 0 0
25 66085 11109 3599 0 0 0 0 0
26 41654 19402 19201 536 0 0 0 0
27 48808 20447 10279 1259 0 0 0 0
28 68983 11810 0 0 0 0 0 0
29 34631 26993 17044 2125 0 0 0 0
30 5561 19065 20682 11364 11597 7963 3301 1260
31 65019 15665 109 0 0 0 0 0
32 19327 26015 29600 5851 0 0 0 0
33 15128 23405 24345 10059 7856 0 0 0
34 76432 4361 0 0 0 0 0 0
35 73829 6080 884 0 0 0 0 0
36 71050 9743 0 0 0 0 0 0
37 38199 22369 15663 4486 76 0 0 0
38 28752 35327 16714 0 0 0 0 0
39 77008 3785 0 0 0 0 0 0
40 49953 27129 3711 0 0 0 0 0
41 79838 955 0 0 0 0 0 0
42 78312 2481 0 0 0 0 0 0
43 32813 23892 17056 6661 371 0 0 0
44 39700 23980 15554 1559 0 0 0 0
45 80793 0 0 0 0 0 0 0
46 48461 27669 4663 0 0 0 0 0
47 12053 36783 14826 15182 1949 0 0 0
48 63151 16980 662 0 0 0 0 0
49 73125 7668 0 0 0 0 0 0
50 79860 933 0 0 0 0 0 0
51 80255 538 0 0 0 0 0 0
52 73946 6847 0 0 0 0 0 0
53 77581 3212 0 0 0 0 0 0
54 76399 4394 0 0 0 0 0 0
55 55217 23951 1625 0 0 0 0 0
56 72564 8229 0 0 0 0 0 0
57 57443 22644 706 0 0 0 0 0
58 78988 1805 0 0 0 0 0 0
59 35315 37942 7536 0 0 0 0 0
60 67881 12912 0 0 0 0 0 0
61 52226 24749 3818 0 0 0 0 0
62 62517 16663 1613 0 0 0 0 0
63 46937 31441 2415 0 0 0 0 0
64 68974 11213 606 0 0 0 0 0
65 80793 0 0 0 0 0 0 0
66 67130 12561 1102 0 0 0 0 0
67 37887 27036 11834 4036 0 0 0 0
68 51139 27047 2205 402 0 0 0 0
69 61791 16492 2510 0 0 0 0 0
70 73231 7354 208 0 0 0 0 0
71 22117 22512 20026 14172 1966 0 0 0
72 10060 26081 33450 7590 3612 0 0 0
73 73696 7097 0 0 0 0 0 0
74 80793 0 0 0 0 0 0 0
75 77528 3008 257 0 0 0 0 0
76 34905 34766 10606 516 0 0 0 0
77 54799 21256 4116 622 0 0 0 0
78 57557 21260 1976 0 0 0 0 0
79 46966 28430 5397 0 0 0 0 0
80 46394 30478 3921 0 0 0 0 0
81 80793 0 0 0 0 0 0 0
82 72153 8640 0 0 0 0 0 0
83 57732 21901 1160 0 0 0 0 0
84 73061 7732 0 0 0 0 0 0
85 78108 2685 0 0 0 0 0 0
86 64841 15952 0 0 0 0 0 0
87 75894 4488 411 0 0 0 0 0
88 77088 3705 0 0 0 0 0 0
89 48328 32465 0 0 0 0 0 0
90 80793 0 0 0 0 0 0 0
91 36647 30720 11100 1868 458 0 0 0
92 54513 24626 1654 0 0 0 0 0
93 73946 6847 0 0 0 0 0 0
94 69274 10406 1113 0 0 0 0 0
95 43213 20293 14643 2644 0 0 0 0
96 65015 15778 0 0 0 0 0 0
97 75428 5365 0 0 0 0 0 0
98 65968 13004 1821 0 0 0 0 0
99 78629 2164 0 0 0 0 0 0
100 52405 15982 9604 2802 0 0 0 0
(For the desktop, the amount of data we generated is too much, so we need to reduce the amount of data we use. To do so we average over time blocks, here of length 100.)
AveragedAbundance <- result$Abundance %>% data.frame(
) %>% dplyr::mutate(
time = floor(time/100)*100
) %>% dplyr::group_by(
time
) %>% dplyr::summarise(
dplyr::across(.fns = ~ mean(.x))
)
We can perform a PCA and see if are data can be summarised by a small number of dimensions. As we shall see, constraining to the first 25 principal components does not harm the system much.
PCA <- prcomp(AveragedAbundance %>% dplyr::select_if(~ any(. > 0)),
center = TRUE, scale. = TRUE, rank. = 25)
There is not actually a lot of dependence within the system it appears. Consider the amount of variance explained by the first six principal components (ordered, as is tradition, by amount of variation explained).
head(summary(PCA)$importance[2, ])
PC1 PC2 PC3 PC4 PC5 PC6
0.03749 0.03273 0.03095 0.02751 0.02501 0.02482
The amount of explained variation is as follows.
sum(summary(PCA)$importance[2, ])
[1] 0.99994
The traditional biplot follows, but despite the seeming presence of patterns, so little of the variance is explained that is probably not worth further examination.
ggplot2::autoplot(PCA, loadings = TRUE,
data = AveragedAbundance %>% dplyr::select_if(~ any(. > 0)),
colour = "time")
We can try again with presence-absence data instead.
AveragedPA <- result$Abundance %>% data.frame(
) %>% dplyr::mutate(
time = floor(time/100)*100
) %>% dplyr::mutate(
dplyr::across(.cols = !time, .fns = ~ .x > 0)
) %>% dplyr::group_by(
time
) %>% dplyr::summarise(
dplyr::across(.fns = ~ mean(.x))
)
PCAPA <- prcomp(AveragedPA %>% dplyr::select_if(~ any(. > 0)),
center = TRUE, scale. = TRUE, rank. = 25)
head(summary(PCAPA)$importance[2, ])
PC1 PC2 PC3 PC4 PC5 PC6
0.05537 0.04584 0.04483 0.03901 0.03578 0.03410
So it is a bit better, but not by much.
ggplot2::autoplot(PCAPA, loadings = TRUE,
data = AveragedPA %>% dplyr::select_if(~ any(. > 0)),
colour = "time")
I should note that this is not a failure of the method persay. What this says to me is that dimension reduction of the system cannot reduce the system to a human-readable set of descriptors. The system is still being reduced (from 100 Species * 10 Environments to 25 Components to cover
sum(summary(PCA)$importance[2, ])
[1] 0.99994
of the variation). My initial hypothesis for when the systems are coupled is that, as coupling increases, the number of principal components should decrease, since one system’s changes will be better predictors of the next system’s changes. (An interesting related question; do the number of principal components correlate with the amount of biodiversity in the system?
ifrm(AveragedAbundance)
ifrm(AveragedPA)
ifrm(PCA)
ifrm(PCAPA)
Since trying to reduce the system dimensionally does not work (well enough) a next attempt might be to consider how the pair-wise beta diversity changes over the course of the simulation. Initial problems include that there are quite a few measures of beta diversity as well as that the (square of the) number of environments will determine how many entries we need to consider.
To try this, we are going to reorganise our data into data frame with attributes of environment, time, and traits/species.
Environments <- lapply(
1:result$NumEnvironments,
function(i, abund, numSpecies) {
time <- abund[, 1]
env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
retval <- data.frame(Environment = toString(i),
Time = time
)
retval <- cbind(retval, env)
colnames(retval) <- c("Environment", "Time",
paste0("Basal", 1:34),
paste0("Consumer", 35:100)
)
return(retval)
},
abund = result$Abundance %>% data.frame(
) %>% dplyr::mutate(
time = floor(time/100)*100
) %>% dplyr::group_by(
time
) %>% dplyr::summarise(
dplyr::across(.fns = ~ mean(.x))
),
numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)
Environments <- dplyr::bind_rows(Environments)
Environments <- Environments %>% dplyr::filter(
dplyr::if_any(dplyr::starts_with(c("Basal", "Consumer")), ~ (.x != 0)),
)
One perhaps interesting idea is to try to group the environments by cluster by considering the abundances (or presence-absence) of the species present as traits. The question then is whether the environments have a tendency to attract to specific points or if they wander around each other without relation. (The latter is more neutral.) If they appear to be convergent (which so far would seem to disagree mostly with the alpha diversity analyses), then that would imply that dynamics determine the majority of the system’s fate, while if they instead wander more randomly, then they would appear to be dominated by the neutral mechanisms. (Of course, this is affected by the rate of the neutral mechanisms, so it exists on a definite continuum.) (Note, of course, that cluster analysis usually includes something like PCA, so we would not necessarily expect a different result.)
My working hypothesis for how to use the distance measures is that we need a metric measure (Kindt and Coe, 2005) with any identification scrubbed off (site, time collected). Kindt and Coe suggest that taking the square root of Bray-Curtis makes it metric and thus usable for mathematical distance analyses (such as hierarchical clustering). The vegan package recomends using Jaccard instead of Bray-Curtis, for the same reason.
EnvironmentDistance <- Environments %>% dplyr::select(
-Environment, -Time
) %>% vegan::vegdist(method = "jaccard")
EnvDistClust <- hclust(EnvironmentDistance)
plot(EnvDistClust, labels = FALSE)
Note that are major breakdowns extremely early yielding a large number of clusters but there are expected to be 10 environments (if historical contingency mattered most) or a small number of clusters (if they are converging). The large number of clusters says to me that the neutral dynamics are jumping rapidly from cluster to cluster.
We can try again with presence absence instead.
EnvironmentDistancePA <- Environments %>% dplyr::select(
-Environment, -Time
) %>% vegan::vegdist(method = "jaccard", binary = TRUE)
EnvDistClustPA <- hclust(EnvironmentDistancePA)
plot(EnvDistClustPA, labels = FALSE)
ifrm(Diversity)
ifrm(EnvDistClust)
ifrm(EnvDistClustPA)
ifrm(EnvDiversity)
ifrm(Environments)
ifrm(EnvironmentDistance)
ifrm(EnvironmentDistancePA)
We repeat the procedures above for the linear set. We use the same pool and interaction matrices, but change the space so that ``adjacent’’ communities can disperse (e.g. 1 <-> 2 <-> 3 <->…<-> 9 <-> 10). Similarly, the history is the same in terms of arrivals and extinctions.
load(file.path(
"..", "experiments", "MNA-FirstAttempt-Result-Env10-Line.RData")
)
toRemove <- result$Abundance <= result$Parameters$EliminationThreshold
result$Abundance[toRemove] <- 0
ifrm(toRemove)
With 10 environments, it is probably not helpful to check 10 individual abundance curves, but looking at the first one might be helpful.
LawMorton1996_PlotAbundance(result$Abundance[seq(from = 1,
to = nrow(result$Abundance),
by = 10), c(1, 2:101)]) -> obj;
obj + ggplot2::scale_y_log10() + ggplot2::guides(color = FALSE)
Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
Warning: Transformation introduced infinite values in continuous y-axis
Every vertical line is a species introduction or extinction by neutral dynamics. Note that, due to dispersal, these should lose some of their abundance to spreading, but if that spreading is faster than the population can recoup its losses, than it can invade and be wiped out (since invadability does not consider spatial dynamics, only local ecological dynamics). The primary notes then are that there are a large number more events appearing to occur (due to the spatial arrangement), but the system as a whole is stabler (since the species that get knocked out can be recaptured from space as well).
Diversity <- lapply(
1:result$NumEnvironments,
function(i, abund, numSpecies) {
time <- abund[, 1]
env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
richness <- rowSums(env != 0)
abundSum <- rowSums(env)
entropy <- env / abundSum
entropy <- - apply(
entropy, MARGIN = 1,
FUN = function(x) {
sum(ifelse(x != 0, x * log(x), 0))
})
species <- apply(
env, MARGIN = 1,
FUN = function(x) {
toString(which(x > 0))
}
)
data.frame(Time = time,
Richness = richness,
Entropy = entropy,
Species = species,
Environment = i)
},
abund = result$Abundance,
numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)
Diversity <- dplyr::bind_rows(Diversity)
Diversity <- Diversity %>% dplyr::mutate(
Evenness = Entropy / log(Richness)
)
ggplot2::ggplot(
Diversity,
ggplot2::aes(
x = Time,
y = Richness,
color = factor(Environment),
alpha = ifelse(Environment == 1, 1, 0.3)
)
) + ggplot2::geom_line(
) + ggplot2::geom_line(
data = Diversity %>% dplyr::group_by(
Time
) %>% dplyr::summarise(
Richness = mean(Richness)
),
mapping = ggplot2::aes(
x = Time,
y = Richness
),
color = "black",
inherit.aes = FALSE
) + ggplot2::guides(
alpha = "none"
) + ggplot2::scale_color_discrete(
"Environment"
)
ggplot2::ggplot(
Diversity,
ggplot2::aes(
x = Time,
y = Entropy,
color = factor(Environment),
alpha = ifelse(Environment == 1, 1, 0.3)
)
) + ggplot2::geom_line(
) + ggplot2::geom_line(
data = Diversity %>% dplyr::group_by(
Time
) %>% dplyr::summarise(
Entropy = mean(Entropy)
),
mapping = ggplot2::aes(
x = Time,
y = Entropy
),
color = "black",
inherit.aes = FALSE
) + ggplot2::guides(
alpha = "none"
) + ggplot2::scale_color_discrete(
"Environment"
)
Warning: Removed 1123 row(s) containing missing values (geom_path).
Warning: Removed 110 row(s) containing missing values (geom_path).
plotly::plot_ly(data = Diversity %>% dplyr::filter(Environment < 2),
x = ~Richness, y = ~Entropy, z = ~Time, type = "scatter3d",
mode = "lines", opacity = 1, line = list(color = ~Time))
ggplot2::ggplot(
Diversity,
ggplot2::aes(
x = Time,
y = Evenness,
color = factor(Environment),
alpha = ifelse(Environment == 1, 1, 0.3)
)
) + ggplot2::geom_line(
) + ggplot2::geom_line(
data = Diversity %>% dplyr::group_by(
Time
) %>% dplyr::summarise(
Evenness = mean(Evenness)
),
mapping = ggplot2::aes(
x = Time,
y = Evenness
),
color = "black",
inherit.aes = FALSE
) + ggplot2::guides(
alpha = "none"
) + ggplot2::scale_color_discrete(
"Environment"
)
Warning: Removed 1440 row(s) containing missing values (geom_path).
Warning: Removed 136 row(s) containing missing values (geom_path).
EnvDiversity <- lapply(
1:((ncol(result$Abundance) - 1) / result$NumEnvironments),
function(i, abund, numSpecies) {
time <- abund[, 1]
env <- abund[, 1 + i + numSpecies * (1:result$NumEnvironments - 1)]
richness <- rowSums(env != 0)
abundSum <- rowSums(env)
environments <- apply(
env, MARGIN = 1,
FUN = function(x) {
toString(which(x > 0))
}
)
data.frame(Time = time,
Richness = richness,
Abundance = abundSum,
Species = i,
Environments = environments)
},
abund = result$Abundance,
numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)
EnvDiversity <- dplyr::bind_rows(EnvDiversity)
ggplot2::ggplot(
EnvDiversity %>% dplyr::filter(Richness > 1),
aes(x = Time, y = Richness, color = Species)
) + geom_point(
alpha = 0.01, size = 3
) + guides(
color = "none"
)
with(EnvDiversity %>% mutate(
Time = floor(Time)
) %>% group_by(
Time, Species
) %>% summarise(
Richness = round(mean(Richness)),
.groups = "drop"
),
table(Species, Richness))
Richness
Species 0 1 2 3 4 5 6 7 8 9 10
1 80793 0 0 0 0 0 0 0 0 0 0
2 80624 0 1 1 2 1 1 2 0 2 159
3 80793 0 0 0 0 0 0 0 0 0 0
4 80047 0 0 2 0 1 0 0 1 0 742
5 79905 0 6 0 6 2 5 1 0 11 857
6 80793 0 0 0 0 0 0 0 0 0 0
7 80793 0 0 0 0 0 0 0 0 0 0
8 80793 0 0 0 0 0 0 0 0 0 0
9 80793 0 0 0 0 0 0 0 0 0 0
10 80793 0 0 0 0 0 0 0 0 0 0
11 80793 0 0 0 0 0 0 0 0 0 0
12 680 0 0 0 0 0 0 0 0 1 80112
13 80639 1 1 1 1 1 1 2 1 1 144
14 84 0 0 0 0 0 0 0 0 0 80709
15 80793 0 0 0 0 0 0 0 0 0 0
16 80793 0 0 0 0 0 0 0 0 0 0
17 80793 0 0 0 0 0 0 0 0 0 0
18 80793 0 0 0 0 0 0 0 0 0 0
19 80793 0 0 0 0 0 0 0 0 0 0
20 80793 0 0 0 0 0 0 0 0 0 0
21 80793 0 0 0 0 0 0 0 0 0 0
22 62 0 0 0 0 0 0 0 0 1 80730
23 80793 0 0 0 0 0 0 0 0 0 0
24 80546 0 1 0 0 0 0 1 0 0 245
25 80793 0 0 0 0 0 0 0 0 0 0
26 80793 0 0 0 0 0 0 0 0 0 0
27 80793 0 0 0 0 0 0 0 0 0 0
28 80793 0 0 0 0 0 0 0 0 0 0
29 80721 1 0 0 1 0 0 0 0 1 69
30 69391 9 2 8 1 8 5 2 8 12 11347
31 80793 0 0 0 0 0 0 0 0 0 0
32 500 0 0 0 0 0 1 0 0 1 80291
33 66491 5 15 10 14 2 14 5 14 17 14206
34 80793 0 0 0 0 0 0 0 0 0 0
35 80332 1 0 1 1 1 1 1 0 2 453
36 80793 0 0 0 0 0 0 0 0 0 0
37 46279 1 1 2 1 1 3 1 2 4 34498
38 80793 0 0 0 0 0 0 0 0 0 0
39 80793 0 0 0 0 0 0 0 0 0 0
40 80793 0 0 0 0 0 0 0 0 0 0
41 80793 0 0 0 0 0 0 0 0 0 0
42 80793 0 0 0 0 0 0 0 0 0 0
43 27630 2 2 0 4 0 0 1 4 7 53143
44 60232 1 8 3 3 7 2 7 5 14 20511
45 80793 0 0 0 0 0 0 0 0 0 0
46 80793 0 0 0 0 0 0 0 0 0 0
47 52191 5 8 3 7 3 7 5 9 12 28543
48 80793 0 0 0 0 0 0 0 0 0 0
49 80793 0 0 0 0 0 0 0 0 0 0
50 80793 0 0 0 0 0 0 0 0 0 0
51 80793 0 0 0 0 0 0 0 0 0 0
52 80793 0 0 0 0 0 0 0 0 0 0
53 80793 0 0 0 0 0 0 0 0 0 0
54 80793 0 0 0 0 0 0 0 0 0 0
55 78735 1 1 2 0 1 2 0 0 2 2049
56 80793 0 0 0 0 0 0 0 0 0 0
57 79020 0 0 1 0 0 0 0 1 0 1771
58 80793 0 0 0 0 0 0 0 0 0 0
59 80793 0 0 0 0 0 0 0 0 0 0
60 80793 0 0 0 0 0 0 0 0 0 0
61 80466 0 1 0 0 1 1 1 0 1 322
62 80793 0 0 0 0 0 0 0 0 0 0
63 80793 0 0 0 0 0 0 0 0 0 0
64 80793 0 0 0 0 0 0 0 0 0 0
65 80793 0 0 0 0 0 0 0 0 0 0
66 80793 0 0 0 0 0 0 0 0 0 0
67 80793 0 0 0 0 0 0 0 0 0 0
68 79023 1 5 4 7 5 6 4 4 4 1730
69 78135 2 2 1 1 1 6 1 2 5 2637
70 80793 0 0 0 0 0 0 0 0 0 0
71 1394 0 0 0 0 0 0 0 0 1 79398
72 292 1 0 0 0 0 1 0 1 0 80498
73 80793 0 0 0 0 0 0 0 0 0 0
74 80793 0 0 0 0 0 0 0 0 0 0
75 80793 0 0 0 0 0 0 0 0 0 0
76 78676 1 4 1 4 0 6 1 3 6 2091
77 79765 0 0 1 0 0 1 0 0 1 1025
78 79253 0 0 1 0 1 0 0 0 0 1538
79 80793 0 0 0 0 0 0 0 0 0 0
80 80793 0 0 0 0 0 0 0 0 0 0
81 80793 0 0 0 0 0 0 0 0 0 0
82 80793 0 0 0 0 0 0 0 0 0 0
83 80793 0 0 0 0 0 0 0 0 0 0
84 80793 0 0 0 0 0 0 0 0 0 0
85 80793 0 0 0 0 0 0 0 0 0 0
86 64862 3 6 4 7 4 6 3 8 6 15884
87 80793 0 0 0 0 0 0 0 0 0 0
88 80793 0 0 0 0 0 0 0 0 0 0
89 80793 0 0 0 0 0 0 0 0 0 0
90 80793 0 0 0 0 0 0 0 0 0 0
[ reached getOption("max.print") -- omitted 10 rows ]
AveragedAbundance <- result$Abundance %>% data.frame(
) %>% dplyr::mutate(
time = floor(time/100)*100
) %>% dplyr::group_by(
time
) %>% dplyr::summarise(
dplyr::across(.fns = ~ mean(.x))
)
PCA <- prcomp(AveragedAbundance %>% dplyr::select_if(~ any(. > 0)),
center = TRUE, scale. = TRUE, rank. = 25)
head(summary(PCA)$importance[2, ])
PC1 PC2 PC3 PC4 PC5 PC6
0.19145 0.11089 0.10844 0.07943 0.05993 0.04530
sum(summary(PCA)$importance[2, ])
[1] 0.99996
ggplot2::autoplot(PCA, loadings = TRUE,
data = AveragedAbundance %>% dplyr::select_if(~ any(. > 0)),
colour = "time")
AveragedPA <- result$Abundance %>% data.frame(
) %>% dplyr::mutate(
time = floor(time/100)*100
) %>% dplyr::mutate(
dplyr::across(.cols = !time, .fns = ~ .x > 0)
) %>% dplyr::group_by(
time
) %>% dplyr::summarise(
dplyr::across(.fns = ~ mean(.x))
)
PCAPA <- prcomp(AveragedPA %>% dplyr::select_if(~ any(. > 0)),
center = TRUE, scale. = TRUE, rank. = 25)
head(summary(PCAPA)$importance[2, ])
PC1 PC2 PC3 PC4 PC5 PC6
0.21750 0.12185 0.08092 0.05943 0.05215 0.04468
ggplot2::autoplot(PCAPA, loadings = TRUE,
data = AveragedPA %>% dplyr::select_if(~ any(. > 0)),
colour = "time")
sum(summary(PCA)$importance[2, ])
[1] 0.99996
ifrm(AveragedAbundance)
ifrm(AveragedPA)
ifrm(PCA)
ifrm(PCAPA)
Environments <- lapply(
1:result$NumEnvironments,
function(i, abund, numSpecies) {
time <- abund[, 1]
env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
retval <- data.frame(Environment = toString(i),
Time = time
)
retval <- cbind(retval, env)
colnames(retval) <- c("Environment", "Time",
paste0("Basal", 1:34),
paste0("Consumer", 35:100)
)
return(retval)
},
abund = result$Abundance %>% data.frame(
) %>% dplyr::mutate(
time = floor(time/100)*100
) %>% dplyr::group_by(
time
) %>% dplyr::summarise(
dplyr::across(.fns = ~ mean(.x))
),
numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)
Environments <- dplyr::bind_rows(Environments)
Environments <- Environments %>% dplyr::filter(
dplyr::if_any(dplyr::starts_with(c("Basal", "Consumer")), ~ (.x != 0)),
)
EnvironmentDistance <- Environments %>% dplyr::select(
-Environment, -Time
) %>% vegan::vegdist(method = "jaccard")
EnvDistClust <- hclust(EnvironmentDistance)
plot(EnvDistClust, labels = FALSE)
EnvironmentDistancePA <- Environments %>% dplyr::select(
-Environment, -Time
) %>% vegan::vegdist(method = "jaccard", binary = TRUE)
EnvDistClustPA <- hclust(EnvironmentDistancePA)
plot(EnvDistClustPA, labels = FALSE)
ifrm(Diversity)
ifrm(EnvDistClust)
ifrm(EnvDistClustPA)
ifrm(EnvDiversity)
ifrm(Environments)
ifrm(EnvironmentDistance)
ifrm(EnvironmentDistancePA)
We repeat the procedures above for the linear set. We use the same pool and interaction matrices, but change the space so that ``adjacent’’ communities can disperse (e.g. 1 <-> 2 <-> 3 <->…<-> 9 <-> 10). Similarly, the history is the same in terms of arrivals and extinctions.
load(file.path(
"..", "experiments", "MNA-FirstAttempt1e+06-Result-Env10-Line.RData")
)
toRemove <- result$Abundance <= result$Parameters$EliminationThreshold
result$Abundance[toRemove[1:(length(toRemove)/2)]] <- 0
result$Abundance[toRemove[(length(toRemove)/2 + 1):length(toRemove)]] <- 0
ifrm(toRemove)
With 10 environments, it is probably not helpful to check 10 individual abundance curves, but looking at the first one might be helpful.
LawMorton1996_PlotAbundance(result$Abundance[seq(from = 1,
to = nrow(result$Abundance),
by = 10), c(1, 2:101)]) -> obj;
obj + ggplot2::scale_y_log10() + ggplot2::guides(color = FALSE)
Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
Warning: Transformation introduced infinite values in continuous y-axis
Diversity <- lapply(
1:result$NumEnvironments,
function(i, abund, numSpecies) {
time <- abund[, 1]
env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
richness <- rowSums(env != 0)
abundSum <- rowSums(env)
entropy <- env / abundSum
entropy <- - apply(
entropy, MARGIN = 1,
FUN = function(x) {
sum(ifelse(x != 0, x * log(x), 0))
})
species <- apply(
env, MARGIN = 1,
FUN = function(x) {
toString(which(x > 0))
}
)
data.frame(Time = time,
Richness = richness,
Entropy = entropy,
Species = species,
Environment = i)
},
abund = result$Abundance,
numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)
Diversity <- dplyr::bind_rows(Diversity)
Diversity <- Diversity %>% dplyr::mutate(
Evenness = Entropy / log(Richness)
)
ggplot2::ggplot(
Diversity,
ggplot2::aes(
x = Time,
y = Richness,
color = factor(Environment),
alpha = ifelse(Environment == 1, 1, 0.3)
)
) + ggplot2::geom_line(
) + ggplot2::geom_line(
data = Diversity %>% dplyr::group_by(
Time
) %>% dplyr::summarise(
Richness = mean(Richness)
),
mapping = ggplot2::aes(
x = Time,
y = Richness
),
color = "black",
inherit.aes = FALSE
) + ggplot2::guides(
alpha = "none"
) + ggplot2::scale_color_discrete(
"Environment"
)
geom_path: Each group consists of only one observation. Do you need to adjust the group
aesthetic?
ggplot2::ggplot(
Diversity,
ggplot2::aes(
x = Time,
y = Entropy,
color = factor(Environment),
alpha = ifelse(Environment == 1, 1, 0.3)
)
) + ggplot2::geom_line(
) + ggplot2::geom_line(
data = Diversity %>% dplyr::group_by(
Time
) %>% dplyr::summarise(
Entropy = mean(Entropy)
),
mapping = ggplot2::aes(
x = Time,
y = Entropy
),
color = "black",
inherit.aes = FALSE
) + ggplot2::guides(
alpha = "none"
) + ggplot2::scale_color_discrete(
"Environment"
)
Warning: Removed 197949 row(s) containing missing values (geom_path).
Warning: Removed 1 row(s) containing missing values (geom_path).
plotly::plot_ly(data = Diversity %>% dplyr::filter(Environment < 2),
x = ~Richness, y = ~Entropy, z = ~Time, type = "scatter3d",
mode = "lines", opacity = 1, line = list(color = ~Time))
ggplot2::ggplot(
Diversity,
ggplot2::aes(
x = Time,
y = Evenness,
color = factor(Environment),
alpha = ifelse(Environment == 1, 1, 0.3)
)
) + ggplot2::geom_line(
) + ggplot2::geom_line(
data = Diversity %>% dplyr::group_by(
Time
) %>% dplyr::summarise(
Evenness = mean(Evenness)
),
mapping = ggplot2::aes(
x = Time,
y = Evenness
),
color = "black",
inherit.aes = FALSE
) + ggplot2::guides(
alpha = "none"
) + ggplot2::scale_color_discrete(
"Environment"
)
Warning: Removed 294140 row(s) containing missing values (geom_path).
Warning: Removed 1 row(s) containing missing values (geom_path).
EnvDiversity <- lapply(
1:((ncol(result$Abundance) - 1) / result$NumEnvironments),
function(i, abund, numSpecies) {
time <- abund[, 1]
env <- abund[, 1 + i + numSpecies * (1:result$NumEnvironments - 1)]
richness <- rowSums(env != 0)
abundSum <- rowSums(env)
environments <- apply(
env, MARGIN = 1,
FUN = function(x) {
toString(which(x > 0))
}
)
data.frame(Time = time,
Richness = richness,
Abundance = abundSum,
Species = i,
Environments = environments)
},
abund = result$Abundance,
numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)
EnvDiversity <- dplyr::bind_rows(EnvDiversity)
ggplot2::ggplot(
EnvDiversity %>% dplyr::filter(Richness > 1),
aes(x = Time, y = Richness, color = Species)
) + geom_point(
alpha = 0.01, size = 3
) + guides(
color = "none"
)
with(EnvDiversity %>% mutate(
Time = floor(Time)
) %>% group_by(
Time, Species
) %>% summarise(
Richness = round(mean(Richness)),
.groups = "drop"
),
table(Species, Richness))
Richness
Species 0 1 4 5
1 1 0 0 0
2 1 0 0 0
3 1 0 0 0
4 1 0 0 0
5 1 0 0 0
6 1 0 0 0
7 1 0 0 0
8 1 0 0 0
9 1 0 0 0
10 1 0 0 0
11 1 0 0 0
12 0 0 0 1
13 1 0 0 0
14 0 0 1 0
15 1 0 0 0
16 1 0 0 0
17 1 0 0 0
18 1 0 0 0
19 1 0 0 0
20 1 0 0 0
21 1 0 0 0
22 0 0 0 1
23 1 0 0 0
24 1 0 0 0
25 1 0 0 0
26 1 0 0 0
27 1 0 0 0
28 1 0 0 0
29 1 0 0 0
30 1 0 0 0
31 1 0 0 0
32 0 1 0 0
33 1 0 0 0
34 1 0 0 0
35 1 0 0 0
36 1 0 0 0
37 1 0 0 0
38 1 0 0 0
39 1 0 0 0
40 1 0 0 0
41 1 0 0 0
42 1 0 0 0
43 1 0 0 0
44 1 0 0 0
45 1 0 0 0
46 1 0 0 0
47 0 1 0 0
48 1 0 0 0
49 1 0 0 0
50 1 0 0 0
51 1 0 0 0
52 1 0 0 0
53 1 0 0 0
54 1 0 0 0
55 1 0 0 0
56 1 0 0 0
57 1 0 0 0
58 1 0 0 0
59 1 0 0 0
60 1 0 0 0
61 1 0 0 0
62 1 0 0 0
63 1 0 0 0
64 1 0 0 0
65 1 0 0 0
66 1 0 0 0
67 1 0 0 0
68 1 0 0 0
69 1 0 0 0
70 1 0 0 0
71 1 0 0 0
72 0 1 0 0
73 1 0 0 0
74 1 0 0 0
75 1 0 0 0
76 1 0 0 0
77 1 0 0 0
78 1 0 0 0
79 1 0 0 0
80 1 0 0 0
81 1 0 0 0
82 1 0 0 0
83 1 0 0 0
84 1 0 0 0
85 1 0 0 0
86 1 0 0 0
87 1 0 0 0
88 1 0 0 0
89 1 0 0 0
90 1 0 0 0
91 1 0 0 0
92 1 0 0 0
93 1 0 0 0
94 1 0 0 0
95 1 0 0 0
96 1 0 0 0
97 1 0 0 0
98 1 0 0 0
99 1 0 0 0
100 1 0 0 0
AveragedAbundance <- result$Abundance %>% data.frame(
) %>% dplyr::mutate(
time = floor(time/100)*100
) %>% dplyr::group_by(
time
) %>% dplyr::summarise(
dplyr::across(.fns = ~ mean(.x))
)
PCA <- prcomp(AveragedAbundance %>% dplyr::select_if(~ any(. > 0)),
center = TRUE, scale. = TRUE, rank. = 25)
Error in prcomp.default(AveragedAbundance %>% dplyr::select_if(~any(. > :
cannot rescale a constant/zero column to unit variance